
# get the command arguments from perl
# rsm number
# runfolder
# run number

args=(commandArgs(TRUE))

for(i in 1:length(args)){
    eval(parse(text=args[[i]]))
}	 

#write.csv(cbind(rsm,runfolder,run),file=paste(runfolder,"/rep",run,".csv",sep=""))



# Plotrix - hex grid plotting
# Ape - phylogenetic tree plotting
library(plotrix)
library(ape)


# overload color2D.matplot
color2D.matplot <- function (x, redrange = c(0, 1), greenrange = c(0, 1), bluerange = c(0, 
    1), extremes = NA, cellcolors = NA, show.legend = FALSE, 
    nslices = 10, xlab = "Column", ylab = "Row", do.hex = FALSE, 
    axes = TRUE, show.values = FALSE, vcol = "white", vcex = 1, 
    border = "black", na.color = NA, ...) 
{
    if (is.matrix(x) || is.data.frame(x)) {
        xdim <- dim(x)
        if (is.data.frame(x)) 
            x <- unlist(x)
        else x <- as.vector(x)
        oldpar <- par("xaxs", "yaxs", "xpd", "mar")
        par(xaxs = "i", yaxs = "i")
        if (do.hex) 
            #par(mar = c(5, 4, 4, 4))
        plot(c(0, xdim[2]), c(0, xdim[1]), xlab = xlab, ylab = ylab, 
            type = "n", axes = FALSE, ...)
        oldpar$usr <- par("usr")
        if (!do.hex) {
            box()
            pos <- 0
        }
        else pos <- -0.3
        if (axes) {
            xticks <- pretty(0:xdim[2])[-1]
            axis(1, at = xticks - 0.5, labels = xticks, pos = pos)
            yticks <- pretty(0:xdim[1])[-1]
            axis(2, at = xdim[1] - yticks + 0.5, yticks)
        }
        if (all(is.na(cellcolors))) 
            cellcolors <- color.scale(x, redrange, greenrange, 
                bluerange, extremes, na.color)
        if (do.hex) {
            par(xpd = TRUE)
            offset <- 0
            if (length(border) < xdim[1] * xdim[2]) 
                border <- rep(border, length.out = xdim[1] * 
                  xdim[2])
            for (row in 1:xdim[1]) {
                for (column in 0:(xdim[2] - 1)) {
                  hexagon(column + offset, xdim[1] - row, col = cellcolors[row + 
                    xdim[1] * column], border = border[row + 
                    xdim[1] * column])
                  if (show.values) 
                    text(column + offset + 0.5, xdim[1] - row + 
                      0.5, x[row + column * xdim[1]], col = vcol, 
                      cex = vcex)
                }
                offset <- ifelse(offset, 0, 0.5)
            }
            par(xpd = FALSE)
        }
        else {
            rect(sort(rep((1:xdim[2]) - 1, xdim[1])), rep(seq(xdim[1] - 
                1, 0, by = -1), xdim[2]), sort(rep(1:xdim[2], 
                xdim[1])), rep(seq(xdim[1], 1, by = -1), xdim[2]), 
                col = cellcolors, border = FALSE)
            if (show.values) 
                text(sort(rep((1:xdim[2]) - 0.5, xdim[1])), rep(seq(xdim[1] - 
                  0.5, 0, by = -1), xdim[2]), round(x, show.values), 
                  col = vcol, cex = vcex)
        }
        naxs <- which(is.na(x))
        xy <- par("usr")
        plot.din <- par("din")
        plot.pin <- par("pin")
        bottom.gap <- (xy[3] - xy[4]) * (plot.din[2] - plot.pin[2])/(2 * 
            plot.pin[2])
        grx1 <- xy[1]
        gry1 <- bottom.gap * 0.95
        grx2 <- xy[1] + (xy[2] - xy[1])/4
        gry2 <- bottom.gap * 0.8
        if (length(cellcolors) > 1) {
            colmat <- col2rgb(c(cellcolors[which.min(x)], cellcolors[which.max(x)]))
            redrange <- colmat[1, ]/255
            greenrange <- colmat[2, ]/255
            bluerange <- colmat[3, ]/255
        }
        rect.col <- color.scale(1:nslices, redrange, greenrange, 
            bluerange)
        if (show.legend) 
            color.legend(grx1, gry1, grx2, gry2, round(range(x, 
                na.rm = TRUE), 2), rect.col = rect.col)
        par(oldpar)
    }
    else cat("x must be a data frame or matrix\n")
}

rsmdata <- read.csv("/home/rumen/cryptsim/analysis/run.simulation.coefficients.models.csv",header=FALSE)
crypts <- read.csv(paste(runfolder,"/rep",run,".crypts",sep=""),header=FALSE)
patterns <- read.csv(paste(runfolder,"/rep",run,".patterns",sep=""),header=FALSE)
glog <- read.csv(paste(runfolder,"/rep",run,".log",sep=""),header=FALSE)


rsmmodel <- rsmdata[rsmdata[,1]==rsm,]
repr <- as.numeric(strsplit(as.character(rsmmodel[1,3])," ")[[1]])
surv <- strsplit(as.character(rsmmodel[2,3])," ")[[1]]
muta <- strsplit(as.character(rsmmodel[3,3])," ")[[1]]

repr.reduced <- repr[repr!=0]
relative.fitness.coefficients <- NULL 
for(i in 1:length(repr.reduced)){
 combos <- combn(repr.reduced, i)
 combos <- combos + 1
 relative.fitness.coefficients <- c(relative.fitness.coefficients, apply(combos, 2, prod))
} 
reproduction.fitness.classes <- c(1.00,sort(unique(relative.fitness.coefficients)))


repro.fitness.multiplier <- t(c(repr[1:5]+1))
repro.fitness <- patterns[,102:106]*repro.fitness.multiplier
repro.fitness[repro.fitness == 0] <- 1
# holds clone id and its reproductive fitness
final.fitness <- cbind(patterns[,1], apply(repro.fitness, 1, prod))


# final.fitness holds: clone_id, fitness, fitness category
final.fitness <- cbind(final.fitness,match(final.fitness[,2],reproduction.fitness.classes))

# Plus one for dead crypts
max.fitness.category <- length(reproduction.fitness.classes)+1

# make wild type be colored with maximum fitness color + 1, and maximum fitness color be colored wild type color - 1
final.fitness[final.fitness[,3]==1,3] <- max.fitness.category

# Color dead crypts with first fitness class
final.fitness <- rbind(final.fitness,c(-1,1,1))

# Get the index of the first and only  occurrence of clone_id from the glog to its matching clone_id from final.fitness table
ind <- match(glog[,2],final.fitness[,1])
glog <- cbind(glog,ind)

# assign the proper colors based on reproductive fitness

mat <- matrix(c(0:89999),nrow=300,ncol=300,byrow=TRUE)

# Get all the unique sampling time points
sampling.time <- unique(glog[,3])

# Set color palette
rgb.palette <- colorRampPalette(c("black", "white"))
palette(rgb.palette(max.fitness.category))

# PNG files to make a movie
dir.create(paste(runfolder,"/images.rep",run,sep=""))
for(i in 1:length(sampling.time)){
png(paste(runfolder,"/images.rep",run,"/rep",run,".",sprintf("%04d",i),".png",sep=""),width=900,height=900,units="px")
par(mar=c(0,0,0,0))
current.time <- sampling.time[i]
current.glog <- glog[glog[,3]==current.time,]
mat.selective.col <- matrix(final.fitness[current.glog[,5],3],nrow=300,ncol=300,byrow=TRUE)
color2D.matplot(mat,cellcolors=mat.selective.col,xlab=NA,ylab=NA,do.hex=TRUE,border=mat.selective.col,axes=FALSE)
dev.off()
}

# mencoder mf://*.png -mf w=900:h=900:fps=2:type=png -ovc copy -oac copy -o output.avi
